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O ; ABSTRACT 

o: 

f^ , For pulsars similar to the one in the Crab Nebula, most of the energy input to 

^ I the surrounding wind nebula occurs on a timescale < 10'^ years; during this time, the 

Qh| nebula expands into freely expanding supernova ejecta. On a timescale ~ 10^ years, the 

O ' interaction of the supernova with the surrounding medium drives a reverse shock front 

-*^ ' toward the center of the remnant, where it crushes the PWN (pulsar wind nebula). One- 

Cd ■ and two-dimensional, two-fluid simulations of the crushing and re-expansion phases of a 

^ . PWN show that (1) these phases are subject to Rayleigh- Taylor instabilities that result 

^ , in the mixing of thermal and nonthermal fluids, and (2) asymmetries in the surrounding 

H I interstellar medium give rise to asymmetries in the position of the PWN relative to the 

pulsar and explosion site. These effects are expected to be observable in the radio 

emission from evolved PWN because of the long lifetimes of radio emitting electrons. 

The scenario can explain the chaotic and asymmetric appearance of the Vela X PWN 

relative to the Vela pulsar without recourse to a directed flow from the vicinity of the 

pulsar. The displacement of the radio nebulae in G327. 1-1.1, MSH15-56 (G326.3-1.8), 

GO.9+0.1, and W44 relative to the X-ray nebulae may be due to this mechanism. On 

timescales much greater than the nebular crushing time, the initial PWN may be mixed 

with thermal gas and become unobservable, so that even the radio emission is dominated 

by recently injected particles. 

Subject headings: ISM: supernova remnants — shock waves — stars: pulsars: general 



1. INTRODUCTION 

Pulsars are expected to be born inside massive stars, so that the evolution of the wind nebulae 
that they produce is expected to depend on a number of different factors: the structure of the 
supernova, the nature of the surrounding medium, the evolution of the pulsar spin-down power, 
and the space velocity of the pulsar. The evolution of pulsar nebulae can be divided into a number of 
phases that are important for their observational appearance (Reynolds & Chevalier 1984, hereafter 
RC84; Chevalier 1998). Initially, the pulsar nebula expands into the freely expanding supernova 
ejecta. The pulsar provides a constant wind power and the swept-up shell of ejecta is accelerated. 
Jun (1998) has carried out two-dimensional simulations of the Rayleigh- Taylor instabilities that 
occur during this phase, with the aim of modeling the Crab Nebula. On a timescale of ~ 10^ yr for 
a pulsar like the Crab, the power input from the pulsar drops steeply so that the nebula expands 
adiabatically. The expansion approaches free expansion within the supernova. 

The next phase of evolution results from the interaction of the supernova remnant with the 
surrounding medium and the inward motion of the reverse shock front that is driven by this in- 
teraction. On a timescale of ~ 10^ yr, the reverse shock makes its way back to the center of the 
remnant. The pulsar nebula that has been created by the early energy injection is compressed 
during this phase. Because of the synchrotron losses of high energy particles, this nebula is best 
observed at radio wavelengths. At the same time, the continued wind power from the pulsar can 
create a nebula, including high energy emission, that is localized to the pulsar. If the pulsar has 
a velocity of lOO's of km s~^, as is quite likely for a normal pulsar (Lyne &: Lorimer 1994), a bow 
shock nebula can form around the pulsar and this nebula can separate from the crushed nebula 
from the earlier phase. 

RC84 made some approximate estimates of the effects of crushing a pulsar nebula by the 
external supernova remnant, with an emphasis on the synchrotron luminosity. Van der Swaluw 
et al. (2000) carried out one-fluid, one-dimensional simulations of the crushing process, finding 
that there are considerable transient effects before the nebula settles into a slow expansion. Our 
aim here is also to investigate the hydrodynamics of the interaction, with attention to instabilities. 
In § 2, we present the basic parameters that are needed to model the pulsar nebula/supernova 
remnant interaction. A difference with the work of van der Swaluw et al. (2000) is that they 
assumed supernova ejecta with a constant density profile, whereas we allow for an outer power law 
profile. Numerical simulations in one and two dimensions are presented in § 3 and § 4, respectively. 
We allow for different adiabatic indices in the relativistic fluid and the thermal gas fluid. In the 
two-dimensional simulations, the effect of a density gradient in the external medium is treated in 
addition to a constant density ambient medium. We believe this situation is relevant to the Vela 
X radio pulsar nebula in the Vela supernova remnant, which is discussed in § 5 along with other 
remnants. The conclusions are in S 6. 



2. A MODEL FOR THE PWN/SNR INTERACTION 

To investigate the dynamical evolution of the pulsar wind nebula/supernova remnant system 
beyond the free expansion phase, we have constructed a simple model based on the expanding 
pulsar bubble solution (Chevalier 1977) and the self-similar driven wave (Chevalier 1982). We 
begin by assuming the stellar material ejected by the Type II supernova is expanding ballistically 
(r = vt) and can be described by a steep outer power law density profile inside of which the density 
is constant. The density profile is parameterized by the total mass of ejecta, M^j, the kinetic energy 
released in the supernova, Egn, and the density power law exponent, n: 

I Av^r^^t^^"^ for r > vd 
""C--**"!^,- for r<nt <" 

where the constant A is given by 

5n — 25 „ c: 
A = -^E..v^^ (2) 

and the velocity at the intersection of the density plateau and the power law is given by 

If we further assume that this ejecta is expanding into a uniform ambient medium with density 
Pa, the expansion of the supernova remnant is described by a self-similar solution given by Chevalier 
(1982). The radius of the forward shock, Ri, in this self-similar driven wave (SSDW) is given by 
(assuming s = in Chevalier's notation) 

Pt^\ ^(n-3)/n (4) 

where the constant a = 1.048 for n = 9, but varies relatively little with the value of n. The 
assumption of a constant density surrounding medium may not generally apply because core collapse 
supernovae have massive star progenitors which are known to affect their surroundings through 
winds and photoionization. The timescale for the reverse shock wave to return to the center (~ 10^ 
years) makes it plausible that the outer shock wave has proceeded to the interstellar medium. If 
a wind bubble has been created by the progenitor star, the inward motion of the reverse shock is 
modified. 

The last component of the model, the pulsar wind nebula, is treated as an adiabatic wind 
bubble driven by a pulsar wind with a kinetic luminosity, Lp. If we use the standard approximation 
of an isobaric bubble bounded by a thin shell of swept up ambient gas, as in Chevalier (1977), 
the expansion of this bubble within the uniform (but decreasing with time) density of the ejecta is 
given by 
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Fig. 1. — The density profile of our model for the interaction of a pulsar nebula with the host 
supernova remnant. The supernova remnant is modeled as a self-similar driven wave (SSDW) 
bounded by a forward shock at Ri and a reverse shock at i?2- The pulsar bubble has swept up a 
thin shell of ejecta at Rp. The edge of the ejecta plateau at Rt is just about to reach R2, after 
which the reverse shock will begin propagating in toward the center. 



This combined model is illustrated in Figure 1 which shows the radial profile of the density 
throughout the structure. The supernova ejecta is separated into three parts: a thin shell at Rp 
that has been swept up by the expanding pulsar bubble, cold freely expanding ejecta between Rp 
and the reverse SNR shock at -R2, and a thin shell of shocked ejecta at R2 that has been decelerated 
by the reverse shock. The discontinuity between R2 and i?i separates the shocked ejecta (on the 
left) and the shocked CSM (on the right). The discontinuity near r = is the termination shock 
of the pulsar wind. 

Several things will happen at an age of ~ 10^ years to change this model. First, the pulsar 
will die out. Without continued power input, the expansion of the pulsar nebula will slow down 
until it becomes frozen into the expanding ejecta. Second, the pulsar nebula will reach the edge of 
the plateau in the SN ejecta and begin accelerating down the steep density gradient. Alternatively, 
the plateau will reach the reverse shock of the SSDW. As the ram pressure at the reverse shock 
rapidly decays away, the reverse shock is driven toward the center of the remnant, resulting in a 
dynamical interaction with the pulsar nebula. 

The exact ordering of these events depends on the parameters of the model described above, 
namely Egn, Mgj, n, pa, Lp and the lifetime of the pulsar. The time of the first event, the "death" 
of the pulsar, is determined by our model for Lp{t). We approximate the power input from the 
pulsar by assuming a constant pulsar magnetic field and braking index p, which yields 

/ ^N -(p+l)/(p-l) 

Lp = Lp^il + -\ , (6) 

where Lpi is the initial pulsar power and p = 3 for the magnetic dipole case. It can be seen that much 
of the power input occurs up to time r, so that r corresponds roughly to the lifetime of the pulsar. 
The initial pulsar power is given in terms of the initial spin period, Pj: Lpi = Iflf/[{p — l)r], where 
/ is the pulsar moment of inertia and Oj = 2tt/ Pi is the initial spin rate. For the Crab pulsar, p is 
measured to be 2.5 and Lp = Ax 10^^ erg s^^, which lead to r = 744 years and Lpi = 3 x 10^^ erg s^^ 
for / = lO'*^ g cm^ (e.g., Chevalier 1977 and references therein). The power input for our model 
C (see Table 1) is considerably above Lpi for the Crab pulsar, but corresponds to Pj ~ 5 ms (for 
p = 3), which is a plausible initial rotation period. 

The second event, the pulsar nebula reaching the edge of the ejecta plateau, Rt = vt*t, occurs 

at 

/ n — 5 \ E^„ , , 

t2 = 2.64 -^ (7) 

\ n J Lp 

assuming constant pulsar luminosity. If, however, the pulsar has already died away, the nebula will 
expand slower and this will be an underestimate to t2- The third event, the inward motion of the 
reverse shock, begins when the edge of the ejecta plateau reaches the reverse shock, i.e., Rt = R2'- 

where the numerical value is quoted for n = 9, for which B = 



Choosing reasonable parameters as listed in Table 2 for Model A, we find ^2 ~ 3, 700 yrs and 
t3 ~ 1, 500 yrs. Thus we expect the reverse shock to begin moving into the interior of the remnant 
where it will encounter the pulsar nebula. For other parameters, however, this sequence may be 
reversed. For example, lowering Egn by 3 would lead to ^2 < ^s; the lower energy means a slower 
expansion, requiring a longer time to sweep up mass but a shorter time for the pulsar to blow 
through the slower ejecta. 



3. ONE-DIMENSIONAL HYDRODYNAMIC SIMULATIONS 

To investigate the dynamical evolution of the pulsar nebula beyond the phase of self-similar 
expansion in a uniform ejecta, we employ numerical simulations that begin shortly before the 
PWN/SNR interaction. These simulations are computed with the VH-1 code, a Lagrangian-Remap 
version of the Piecewise Parabolic Method. The standard VH-1 code was adapted to treat two fluids 
with two different values of 7 such that the pulsar wind can be modeled with 7 = 4/3 while the 
supernova ejecta and circumstellar gas are modeled with 7 = 5/3. The one-dimensional PWN/SNR 
model is evolved on a grid of 1000 radial zones that expands to follow the forward shock. Several 
simulations on a grids of 750 and 1500 zones produced almost identical results. 

The initial conditions for the SSDW are taken from the solution described in Chevalier (1982), 
scaled to the parameters listed in Table 2 for Model A. The initial conditions for the pulsar nebula 
are taken from a separate hydrodynamic simulation of a wind blown bubble in which the wind is an 
adiabatic gas with 7 = 4/3, a mass loss rate M^ = 2Lp/vi^ and a wind velocity, Up, scaled to yield 
a sound speed in the shocked gas of c/\/3, appropriate for a relativistic gas. The pulsar luminosity 
decays away according to eq. (6) assuming a braking index of p = 3. 

A more realistic model for the pulsar wind nebula would allow for a relativistic, MHD (mag- 
netohydrodynamic) flow, as in the steady state model for the Crab Nebula by Kennel & Coroniti 
(1984a, b). However, once the relativistic wind passes through the wind termination shock, the flow 
velocity drops to c/3. The bulk velocity is not relativistic in the nebula, although the shocked par- 
ticles are highly relativistic. This situation can be approximated by a 7 = 4/3 fluid. The magnetic 
field may be a more significant omission because the magnetic field in the Crab Nebula is thought 
to be in approximate equipartition with the particles. In particular, the magnetic field could affect 
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the growth of instabihties at the edge of the pulsar nebula (Jun, Norman, & Stone 1995; Hester 
et al. 1996). The tangential magnetic field may be especially strong at the edge of the nebula 
(Emmering & Chevalier 1987). The magnetic field may inhibit the growth of the instability along 
the direction of the field, but has less of an effect across the field. 
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Fig. 2. — Evolution of the pulsar nebula within the host supernova remnant as computed with 
the one-dimensional hydrodynamics model. Both the density (top) and the pressure (bottom) are 
depicted in logarithmic scales, with yellow corresponding to high values and purple to low values. 
The radial extent of the grid is scaled to follow the forward shock front, such that the location of 
the forward shock is a horizontal line near the top of each plot. 

The evolution of the PWN/SNR model is illustrated with a space-time diagram in Figure 2. 
The evolution begins at t ~ 2t, just prior to the end of the SSDW phase at ^3. One can identify 



two bands of high-density gas in the early evolution: the pulsar nebula shell made up of swept 
up ejecta, and a shell of shocked ejecta on the trailing side of the SSDW. The simulation begins 
with the edge of the ejecta plateau reaching the reverse shock. The reverse shock and the shell of 
shocked ejecta then begin to push inward, but the SSDW does not have much time to react to the 
dropping ejecta density before the pulsar shell strikes it. 

At an age of about ~ 3000 years the pulsar nebula shell collides with the dense SSDW shell. 
From this point on all of the ejecta is contained in a relatively thin shell (width ~ 0.1 i?i) bounded 
by low density shocked pulsar wind on the inside and by shocked circumstellar gas on the outside. 
This initial impact generates a strong shockwave transmitted through the SSDW. When this shock 
catches up with the forward shock there is an abrupt, but brief, acceleration of the forward shock, 
followed by a faster deceleration than in the SSDW phase. A weak reflected shock is also generated, 
which bounces back and forth through the pulsar bubble at the relatively high sound speed of c/ v3. 
This bouncing wave can be seen as nearly vertical lines in the pressure image of Figure 2. 

Following the collision of the pulsar and reverse SNR shocks, the thin shell of ejecta is gradually 
pushed inward with respect to the outer shockwave by the high pressure of the shocked CSM. 
Initially the pressure in the pulsar bubble is negligible, and the presence of the PWN has no 
effect on the dynamics of the crashing reverse shock. But as the pulsar bubble is compressed, 
the interior pressure rises. For our nominal PWN/SNR parameters, the pressure in the bubble 
becomes comparable to the pressure behind the reverse shock once the PWN has been crushed 
to approximately half the radius of the SNR. Eventually the PWN reaches a minimum radius of 
0.25 -Ri at an age of ~ 35,000 years. At this point the pressure of the compressed pulsar bubble 
accelerates the shell of shocked ejecta back out to a radius of 0.3 Ri. From this point on all of the 
ejecta remains in a thin shell with a width of ~ 10% of Ri, slowly oscillating about an average 
radius of ~ 0.27 i?i. The simulations of van der Swaluw et al. (2000) also show oscillations after 
the crushing phase. 

The extent to which the pulsar bubble is compressed depends on the pressure in the bubble 
compared to the pressure behind the reverse shock as it propagates in toward the pulsar nebula. 
Assuming a constant pulsar luminosity and using the pressure behind the reverse shock in the 
SSDW phase, we find a ratio of pressures at a time of ^3: 



5 = 0.166 l2/5 M3f «-2/i5 R-3/5, 



assuming an ejecta exponent of n = 9. This ratio is ~ 0.2 for our nominal parameters. Note that 
both of these pressures will be slightly lower: Pf, because the pulsar has likely died out by the time 
of the collision, P2 because the pressure in the driven wave will drop once the reverse shock starts 
falling to the center. 

We have run additional models in order to show the effects of changing the energy in the pulsar 
bubble on the evolution of the PWN/SNR interaction. In Model B we dropped the pressure ratio 
in eq. (9) by a factor of '^ 3 through a smaller ejecta mass (4 Mq) and smaller pulsar luminosity 
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Fig. 3. — Evolution of the pulsar nebula through the return of the reverse SNR shock for the three 
different models listed in Table 1: Model A (solid), Model B (long dash), and Model C (short dash). 
The radius of the pulsar nebula is scaled to the (increasing) radius of the SNR as in the previous 
figure. 
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(10^^ erg s^^). In Model C the pressure ratio was increased by a factor of 2.4 by doubling the 
ejecta mass (16 Mq) and increasing the pulsar luminosity (5 x lO^'' erg s~^). The evolution of all 
three models is shown in Figure 3. As expected, the lower the pressure in the pulsar bubble, the 
more the bubble is compressed by the reverse shock. In the lower pressure model the radius of the 
bubble after collapse is only ~ 0.06 i?i. 

The evolution in the high-pressure model is slightly different since these parameters switch 
the ordering of ^2 and t^. In this case the pulsar nebula hits the edge of the ejecta plateau at 
~ r, well before the PWN/SNR collision. There is thus an intermediate phase when the pulsar 
shell accelerates down the steep ejecta gradient before impacting the SSDW. This results in a more 
violent collision, which compresses the SSDW by a factor of ~ 2 and accelerates the forward shock 
considerably more than in the standard model. Despite this more violent beginning, however, the 
subsequent evolution is very similar to the standard model, albeit leaving behind a larger pulsar 
bubble. 

It can be seen in Figure 3 that after the PWN has been crushed, the ratio Rp/Ri evolves 
slowly with time, if the pressure wave induced oscillations are averaged out. If approximate pressure 
equilibrium is established between the PWN and the interior of the blast wave, we expect Rp oc t^'"^ 
(RC84) . The radius of a Sedov blast wave evolves as 

so the slow evolution of Rp/ Ri is expected. Van der Swaluw & Wu (2001) have emphasized this 
point, noting that equating the pressure in the PWN with the interior pressure in the blast wave 
leads to 

^ = im{^^\'\ (11) 

where Epyjn{t) is the internal energy in the PWN at time t. We now assume that the pulsar power 
can be assumed to be constant up to time tp, after which it falls to 0. The value of Epwn{t) can 
be related to the total spindown energy of the pulsar Egd = Lpitp for p = 3. During the initial 
expansion phase {t < tp), the radius of the PWN is given by eq. (5) and Ep^nit) = 5Lpt/ll (RC84). 
For t > tp, the evolution of the PWN is adiabatic, so that Ep^nit) = Ep^nitp) Rp{tp) / Rp{t) , where 
Epwn{tp) = 5£'s(i/ll. With the use of eqs. (5) and (10), eq. (11) can be expressed as 
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Application of eq. (12) to the 3 models discussed here yields Rp/Ri = 0.20 (Model A), 0.11 (Model 
B), and 0.29 (Model C) at t = 5 x 10^^ s if we take tp = r. It can be seen from Fig. (3) that these 
values are slightly lower than seen in the simulations. This discrepancy is due in large part to the 
time-dependence of Lp, such that some pulsar energy continues to be added to the bubble after a 
time r. In simulations with a constant value of Lp up until tp and zero thereafter, the radius of the 
crushed PWN at late times was better fit by this approximation. 
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We have Egd = I^'i/'^, so that the observed value of Rp/Ri can be related to the initial spin 
rate of the pulsar (van der Swaluw & Wu 2001); since Rp/Ri oc ^f^, the value of Rp/Ri appears to 
be more sensitive to Qi than to the other uncertain parameters. However, there is some sensitivity 
to the power input timescale, r, because if r is small, the wind power is injected at an early time 
and the bubble energy is more subject to adiabatic expansion losses. In the magnetic dipole model 
with braking index p = 3, r in eq. (6) is given by (e.g., Shapiro & Teukolsky 1983) 

(13) 



^fBl.Rl^ sin2 a 



where Bns is the magnetic dipole field strength of the pulsar, Rns is its radius, and a is the angle 
between the dipole field and the rotation axis of the pulsar. Ordinary radio pulsars appear to 
have magnetic field strengths in a fairly narrow range around 3 x 10"*^^ G (e.g., Stollman 1987), 
comparable to the field strength of the Crab pulsar. There may be a separate class of strongly 
magnetized neutron stars, magnetars, with Bns ^ 10^"^ G (Duncan & Thompson 1992); these 
objects are expected to rapidly deposit their rotational energy into the surroundings. Turning 
to the initial spin rate and taking r oc O^ as given by eq. (13), we have Rp/Ri oc O^'^. The 
insensitivity of Rp/Ri to Vti is due to the fact that a pulsar with a large spin down energy tends 
to deposit the energy into the nebula early, when it is subject to adiabatic expansion losses. In 
addition, early radiative losses can decrease the energy in the pulsar bubble (see also van der Swaluw 
& Wu 2001). The radiative luminosity of the Crab Nebula at present is ~ 1/3 of the spin-down 
power (Davidson & Fesen 1985), and the higher magnetic field at earlier times can make the losses 
even more important. Atoyan (1999) has described a model for the Crab Nebula in which the 
initial period is ~ 4 msec and most of the spin down energy is lost to radiation. For these combined 
reasons, the observed value of Rp/Ri for an evolved pulsar wind nebula is unlikely to lead to a 
reliable estimate for the initial spin rate. 

To illustrate these effects we show, in Figure 4, the evolution of Rp/Ri for three models with 
the same SNR parameters, but different pulsar winds. In particular, we include two models with 
the same total spin down energy: Model D with Lp = 5x W^^ erg s~^ and Tp = 500 yrs, and Model 
E with Lp = W^^ erg s~^ and Tp = 2500 yrs. Eq. (12) predicts that the pulsar with a more gradual 
power output will produce a PWN at late times that is 50% larger than a PWN from a pulsar with 
the same spin down energy, but releasing it away 5 times faster. 

In addition to the PWN, radiative cooling can be important for the outer SNR, resulting in 
the formation of a dense shell (Chevalier 1974; Cioffi, McKee, & Bertschinger 1988). The shell 
radius expands approximately as Ri oc t^'"^. If both the PWN and the hot interior SNR gas evolve 
adiabatically and R^ <^ Rf, pressure equilibrium between the PWN and the SNR implies that 
Rp oc R^ or Rp/Ri oc t^-^"^^. The size of the PWN nebula relative to the SNR again remains 
approximately constant with time, although now the PWN expands slightly more rapidly. The hot 
SNR gas may not evolve adiabatically but may be subject to radiative cooling just inside the dense 
shell (Chevalier 1974). This has the effect of further increasing the ratio Rp/Ri. 
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Fig. 4. — Evolution of the PWN for models A, D, and E listed in Table 2. The radius of the pulsar 
nebula is scaled to the (increasing) radius of the SNR as in the previous figure. 



4. TWO-DIMENSIONAL HYDRODYNAMIC SIMULATIONS 

4.1. Evolution in a Uniform Medium 

To investigate the effects of instabilities on the morphology of the pulsar nebula, we have 
extended our models to two dimensions. We computed these models on a cylindrical (R-z) grid 
of 2000 X 2000 zones. We assume reflection symmetry about the equator and compute only one 
hemisphere. The initial conditions are taken from the first stage of the one-dimensional simulations 
just after the drop of the pulsar wind termination shock. In order to make these simulations more 
computationally efficient, we have raised the density in the pulsar bubble by a factor of 30. This 
artificially lowers the sound speed in the hot bubble allowing the use of a larger time step. Although 
this increases the time needed to equilibrate the pressure in the bubble, the sound crossing time 
is still relatively short compared to the expansion time, and one-dimensional simulations with and 
without this modification produced identical global evolutions. Again, we evolve this model on a 
moving grid that follows the outer shock front, so that the full resolution of the grid is always used. 

The outer shock front remains very nearly spherical and follows the time evolution of the one 
dimensional model to within 0.1%. The key difference is the instability of the thin shell of ejecta and 
the rapid mixing driven by this instability. Even before the PWN/SNR interaction, the instability 
of the SSDW (Chevalier et al. 1992) begins to spread out the shell of shocked ejecta. However, given 
the short time between the beginning of the simulation and the crash of the reverse shock, there 
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Fig. 5. — Crushing of the pulsar nebula in a two-dimensional simulation using the parameters for 
Model A listed in Table 2. The images show the gas density on a logarithmic scale covering 2.5 
orders of magnitude. The length scale of the images are normalized to take out the expansion of 
the outer shock front. 



is not enough time for this instability to grow to significant amplitude. Following the PWN/SNR 
collision, the deceleration of the ejecta shell by the shocked ambient medium (responsible for the 
SSDW instability) increases in magnitude, thereby driving a more rapid Rayleigh- Taylor instability. 
The result is a broadened shell of mixed ejecta and ambient medium as seen in the first two frames 
of Figure 5. 

Much more dramatic, however, is the Rayleigh- Taylor instability working in the opposite di- 
rection when the high pressure of the compressed pulsar bubble begins to accelerate the ejecta shell 
back outwards. From the ID simulation shown in Figure 2, we see that once Rp/Ri drops below 
~ 0.5, the pressure in the compressed PWN exceeds that in the SNR and the PWN begins to 
decelerate the shell of shocked ejecta. This acceleration of the dense ejecta gas by the low-density 
PWN gas is subject to the Rayleigh- Taylor instability. As a result of this instability, much of the 
ejecta gas continues toward the center of the SNR almost unabated. In the absence of a PWN, 
a spherical reverse shock would reach the center of the remnant at an age of ~ 35, 000 yrs, and 
indeed we see in Figure 5 that some has reached the center by this time. In the same process some 
relativistic gas is displaced from the center. 

We see two competing effects of this instability in the PWN crushing phase. First, as relativistic 
gas is displaced from the center, it escapes the full compression of the infalling ejecta seen in ID. 
Second, the vigorous turbulence driven by this instability leads to rapid mixing of the thermal and 
relativistic gasses. To illustrate these effects, we compare the effective value of Rp/Ri in the ID and 
2D simulations in Figure 6. We compute the radius of the pulsar bubble in the 2D simulation by 
summing up the volume of gas with 7 < 1.66 (i.e., including the partially mixed gas) and calculate 
an effective radius assuming a spherical volume. Prior to the bounce, the ID and 2D simulations 
are nearly identical. However, at the moment of bounce, the PWN in the 2D simulation has a 
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volume twice that of the PWN in the ID simulation. The compressed PWN quickly rebounds in 
the ID simulation, but in 2D the volume of relativistic gas continues to shrink due to numerical 
mixing. 
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Fig. 6. — Evolution of the pulsar nebula through the crash of the reverse SNR shock in one (solid 
line) and two (dashed lines) dimensions with different numerical resolutions. 

This rapid depletion of the volume of the PWN in the 2D simulation is an artifact of numerical 
diffusion; when one numerical zone contains both relativistic gas (7 = 4/3) and ejecta gas (7 = 5/3), 
the mass- weighted average 7 is dominated by the high density of the ejecta. As a result, the mixing 
- and subsequent loss of PWN volume - is strongly dependent on the numerical resolution of the 
simulation. Higher spatial resolution leads to less numerical mixing across the contact interface 
between the two gasses, but it also leads to more vigorous hydrodynamic mixing by resolving 
the turbulence generated on smaller scales. The result of these competing effects can be seen 
in Figure 7, where we compare the distribution of relativistic gas in three 2D simulations with 
varying resolutions. We also plot the effective value of Rp/Ri for these three simulations in Figure 
6. Up until a time of ~ 35, 000 yrs, all three simulations are quite similar. Once the Rayleigh- 
Taylor instability kicks in, however, the resolution-dependent effects of mixing become quite evident. 
Independent of the numerical resolution, it is clear that the mixing of relativisitc gas and thermal 
gas is very efficient. 
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Fig. 7. — The rapid mixing of the thermal gas of the SN ejecta (white) and the relativistic gas of 
the PWN (black) is seen in these images of 7 at a time of 50,000 yrs, shortly after the bounce of 
the pulsar bubble in the corresponding one-dimensional simulation. The three images correspond 
to increasing numerical resolution, with the number of zones per dimension, N = 500 (left), 1000 
(center), 2000 (right). 

4.2. Evolution in Nonuniform Media 



We have repeated our 2D simulation with the addition of a density gradient in the ambient 
medium with the goal of understanding the displacement of the pulsar bubble seen in the Vela 
SNR. Following Dohm-Palmer & Jones (1996), we use a smooth density distribution in the vertical 
direction given by 



P{z) = Pc 



1 + 



■ tanhiz/H] 



(14) 



such that the density contrast, from a maximum at large positive z to a minimum at large negative 
z, is given by (x — 1)^^, and H is the characteristic length scale over which the density changes. 
The simulations shown here use x = 1.2 corresponding to a density contrast of 5. These simulations 
require the computation of the full domain in z (i.e., no assumption of equatorial symmetry). As a 
result, our standard grid of 2000 zones represents only half the resolution used in the previous 2D 
model. 

In Figure 8 we show the evolution of the PWN/SNR system for an ambient medium length 
scale oi H = 1 X 10^^ cm, which corresponds to the size of the SNR at an age of ~ 500 yrs. Early 
in the evolution, we can estimate the asymmetry in the SNR by assuming only radial motion, such 
that the local radius is determined by eq. (4) with the local preshock density. The full density 
contrast of 5 thus produces a relatively minor ratio in radii at the poles of 5^'^ ~ 1.2. The ratio 
of pressures behind the forward shocks is more pronounced: 5^'^ ~ 3.5 Thus, while the forward 
shock is only mildly asymmetric, the interior pressure variation drives a substantial asymmetry in 
the process of crushing the PWN. 

If the ambient medium had an exponential density distribution instead of a step in density. 
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20,000 yrs 



6,000 yrs 



Fig. 8. — The evolution of the pulsar nebula/supernova remnant for the model with a scale height 
of H = 1 X 10^^ cm in the ambient medium. The logarithmic density scale covers 2.5 orders of 
magnitude. 

the distortion of the outer shock front would not be as apparent (e.g., Hnatyk &: Petruk 1999). In 
Figure 9 we show the results of two additional simulations with larger values of H, but all at the 
same age of ~50,000 yrs. In the simulation with H = 1 x 10^'^ cm, the asymmetry in the CSM 
at the time of the PWN/SNR interaction is less than a factor of 2, and the crushed PWN is only 
slightly off center. But a slightly steeper density gradient, H = 3 x 10^^ cm, begins to affect the 
SNR prior to the PWN/SNR interaction, and the crushed PWN is displaced from the center of the 
remnant by roughly 40% the radius of the SNR. 




Fig. 9. — The remnants of the pulsar nebula as defined by the adiabatic index (black is gas with 
7 < 5/3) for simulations with a scale height oi H = Ix 10^^ cm (left) and H = 3x 10^^ cm (center) 
and H = 10^^ cm (right). The solid contours mark the location of the shock front. These snapshots 
are from a time of ~ 50, 000 years, when the SNR has a radius of ^ 1.5 x lO^'' cm. 
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5. COMPARISON WITH OBSERVATIONS 

5.1. The Vela Remnant 

A pulsar nebula that is likely to be in the post-reverse shock phase is Vela X in the Vela 
supernova remnant. This nebula was first identified as a region with a flat radio spectral index 
near the Vela pulsar (Weiler & Panagia 1980). More recent radio studies have shown that some 
of the emission is in filamentary structures as well as a diffuse background (Milne 1995; Frail et 
al. 1997; Bock, Turtle, & Green 1998). The extended Vela X emission is largely to the south of 
the pulsar. The radio observations of Bock et al. (1998) give an especially detailed view of the 
Vela X nebula. The overall appearance is chaotic, in keeping with the expected instabilities during 
the crushing and re-expansion phases. There are a number of radio filaments, and the magnetic 
field lines lie along the filaments; this is expected if the Rayleigh- Taylor instability occurs primarily 
across the magnetic field rather than along it. In their Fig. 5, Bock et al. (1998) show the projected 
velocity of the Vela pulsar, including the pulsar position 12,000 years ago; it is near the North edge 
of the Vela X nebula. Even if this age is incorrect, the direction of the pulsar velocity cannot take 
it to the center of the nebula. We propose that the reason for the asymmetry is an asymmetry in 
the medium surrounding the supernova remnant. 

The interstellar medium surrounding the Vela remnant is known to be inhomogeneous. From 
X-ray observations, Bocchino, Maggio, & Sciortino (1999) deduced that the presence of various 
interstellar components, the lowest density component having a hydrogen density nn ~ 0.06 cm~^. 
From HI observations, Dubner et al. (1998) deduced a preshock density nu ~ 1 cm~^; this denser 
gas is primarily on the north side of the Vela remnant. We thus propose that the supernova occurred 
at the pulsar position (corrected for its velocity), that the higher density to North caused greater 
SNR expansion to the south, and that the asymmetry in the reverse shock front pushed the pulsar 
nebula to the south of the pulsar. The age of the remnant can be estimated from the spindown age 
of the pulsar, tgd = P/[{p — ^)P]i where p is the pulsar braking index, P is the pulsar period, and 
P is its time derivative. The spindown age is the actual age if p is constant, the pulsar magnetic 
field is constant, and the current spin rate is much less than the initial rate. If this last assumption 
is violated, the age is < tsd- For the magnetic dipole model with p = 3, we have tg^ = 11,000 
years for the Vela pulsar (Reichley, Downs, & Morris 1970). However, Lyne et al. (1996) found 
evidence for p = 1.4 it 0.2, which leads to a larger value of tgd ~ 56,000 years. These estimates 
show consistency with the hypothesis that the reverse shock front has recently crushed the pulsar 
nebula in Vela. Another expectation of the reverse shock model is that the pressure in the pulsar 
nebula should be comparable to the pressure in the hot gas in the Vela remnant. Estimates of the 
pressure in the relativistic fiuid from minimum energy synchrotron emission arguments show that 
this is approximately the case (Frail et al. 1997). 

With the standard assumptions, the evolution of the pulsar spin rate is 

/ ^ \ -i/(p-i) 
VL = VLi[l + --\ , (15) 
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where rij is the initial pulsar spin rate and r is approximately the tiniescale for power input by 
the pulsar (see § 3). The present period of the Vela pulsar is P = 2ti/Vl = 89 ms. If p = 1.4 and 
t >> r, we have O k, Q,i{t / t)~'^'^ and the rapid evolution of the spin rate requires a value of the 
initial spin rate that is faster than plausible if r ~ 500 years. This is true even for our model E 
with T ~ 2500 years. However, the value of p = 1.4 by Lyne et al. (1996) is uncertain because of 
the frequent glitching of the pulsar and a value closer to the magnetic dipole value of 3 is possible. 
Of the models presented in § 3, model A appears most suitable because the radius of the pulsar 
nebula is 1/4 of the outer shock radius after the crushing has occurred; the diameter of Vela X is 
~ 2° (Bock et al. 1998), as compared to the ~ 8° diameter of the Vela SNR. 

The fact that the Vela pulsar is to one side of Vela X has led to previous discussions of how 
power might flow from the pulsar into the Vela X region. Observations with ROSAThaYe shown an 
X-ray filament to the south of the pulsar that has been interpreted as emission from a jet extending 
from the pulsar (Markwardt & Ogelman 1995). It is not yet clear whether the X-ray emission is 
thermal or nonthermal (Markwardt & Ogelman 1997). Frail et al. (1997) suggested that some of 
the radio structures are related to the X-ray "jet." However, the radio filament near the X-ray 
feature appears to be similar to other radio filaments in the Vela X region (Bock et al. 1998). Also, 
the X-ray spectrum of the feature is similar to that of the more diffuse X-ray emission in the Vela 
X region (Markwardt & Ogelman 1997). Wijers & Sigurdsson (1997) noted difficulties with the 
high momentum in the proposed X-ray jet and suggested a different model combining fallback onto 
the neutron star with a trail left by a high velocity binary companion. 

The model described here suggests a different interpretation of the "jet" features. As mentioned 
above, the linear feature in the radio synchrotron emission can be the result of the Rayleigh- Taylor 
instability during the crushing process. The instability also results in hot, thermal gas being 
mixed into the region initially occupied by the relativistic fluid. Because of the sweeping action of 
the pulsar nebula, the gas immediately outside the pulsar nebula has a higher density than that 
occupying most of the SNR (Figs. 1 and 2). It is thus plausible that some of the gas mixed into the 
pulsar nebula has a higher emissivity and appears as an X-ray filament. Frail et al. (1997) found 
that the radio and X-ray emission are not coincident but lie near each other, as would be expected 
in our model. 



5.2. Other Remnants 

The return of a reverse shock front to the center of a supernova remnant is not generally 
expected to be spherically symmetric because of inhomogeneities in the surrounding medium. Thus 
we might expect our models to be broadly applicable to composite supernova remnants which have 
both an apparent pulsar nebula and a well-developed shell (see Helfand 1998 for a list of objects). 
The radio emitting electrons are long lived and are the best tracer of the power input from the 
pulsar during its phase of high power output. The X-ray emitting electrons are shorter lived and 
provide an indication of the current position of the pulsar, even if a pulsar has not been directly 
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observed. However, an observable shell may be present before the reverse shock crushing phase. 
This appears to be the case, for example, in the Large Magellanic Cloud remnant 0540-69, which 
has a pulsar nebula very similar to the Crab Nebula (Kaaret et al. 2001). At late times, the mixing 
of the radio emitting electrons with the thermal gas may lead to the disappearance of the radio 
nebula from the initial injection phase. In this case, only emission from near the pulsar, likely to 
be in a trail because of the pulsar velocity, is present. 

In the Vela remnant, the pulsar is surrounded by an X-ray PWN, while the radio emitting PWN 
is displaced from the X-ray emission. A similar situation appears to be present in the remnants 
MSH15-56 (G326.3-1.8; Plucinsky 1998; Gaensler 2001), G327.1-1.1 (Sun, Wang, & Chen 1999), 
and GO. 9+0.1 (Gaensler, Pivovaroff, & Garmire 2001) although pulsars associated with the X-ray 
emission have not yet been discovered in these cases. The broad regions of radio emission are 
excellent candidates for wind nebulae injected early in the pulsar evolution and crushed by the 
reverse shock. 

A less clearcut case is the PWN inside W44 (Frail et al. 1996; Harrus et al. 1996). The radio 
nebula has a bow shock appearance, but the trailing emission is considerably broader and stronger 
than that near the pulsar. The ratio of PWN radius to outer shock radius is ~ 0.1. Although 
smaller than the case of Vela, this situation is expected in our models (see Fig. 3). The age of W44 
is ~ 20, 000 — 30, 000 years, but the relatively high density of its surroundings places it in a more 
evolved category than Vela (Chevalier 1999). Another remnant with relatively dense surroundings 
is IC 443, with an estimated age of 30,000 years (Chevalier 1999). Early radio spectral observations 
indicated the presence of a PWN (Green 1986), and recent X-ray and radio observations clearly 
show a bow shock nebula with a compact X-ray source at the head (Olbert et al. 2001). In this 
case, the radio emission region is narrow and is strongest close to the compact X-ray source; there is 
no clear sign of emission from a remnant nebula that has been crushed by the reverse shock front. 
The trailing emission is an approximately linear feature, but does not point toward the center 
of the remnant. This appears to be another case of an off-center explosion due to asymmetries 
in the surrounding medium. This situation may account for some of the high velocities deduced 
for pulsars from their positions relative to what appears to be their parent SNR (Frail, Goss, & 
Whiteoak 1994). 



6. DISCUSSION AND CONCLUSIONS 

We have considered the evolution of PWN (pulsar wind nebulae) inside of supernova remnants, 
with special attention to the effect of the reverse shock front on the PWN. It is possible that 
synchrotron losses are important for the PWN at early times, but they are unlikely to be significant 
at the times of interest here and the effect can be approximately accounted for by renormalizing 
the pulsar spin down energy. We do not treat the instabilities that can occur early in the evolution 
as a result of the acceleration of the shell outside the pulsar nebula (Jun 1998); these would only 
further contribute to the instabilities that we find later in the evolution. 
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Our simulations allow for two fluids: a 7 = 4/3 fluid for the PWN and a 7 = 5/3 fluid for 
the surrounding gas. We do not allow for the magnetic field that is expected to be mixed in 
with the relativistic particles and that can influence the evolution through magnetic tension. The 
field might affect the Rayleigh- Taylor instabilities that occur when the PWN is compressed by the 
reverse shock gas and subsequently re-expands. The magnetic field can dampen the instability along 
magnetic field lines, but not across them, leading to the creation of filaments along the magnetic 
field direction. There is mixing between the fluids so that the separation of the two fluids depends 
on the numerical resolution of the two-dimensional simulations. The degree of mixing may also 
depend on the magnetic field. 

In addition to causing mixing of thermal gas with the PWN, the reverse shock wave is unlikely 
to be spherically symmetric and can thus push the PWN off center from the explosion site. The 
position of a pulsar can be displaced from the explosion site both because of this effect and because 
of the velocity that the pulsar acquires at birth. 

In comparing the model results with observations, it is important to correctly identify the 
evolutionary state of a PWN. During the initial phase when the PWN is interacting with the freely 
expanding supernova ejecta, the pulsar is expected to be approximately centrally located within 
the PWN. Even if the pulsar has a significant velocity, it will not have a particular position in the 
uniformly expanding medium unless there is a large density gradient in the medium. Nebulae likely 
to be in this phase include the Crab Nebula and 3C 58. 

The next event in the evolution is the cessation of power input from the central pulsar. The 
outer edge of the nebula continues its free expansion, leading to a low surface brightness nebula 
that can have a substantial filling fraction within the external shock wave. The nebula MSH 15-52 
around the pulsar PSR B1509-58 (Seward et al. 1984) may be in this phase. 

Here, we have emphasized the late evolutionary phase after the reverse shock front comes back 
to the center of the remnant. An excellent candidate for this phase of evolution is the Vela X radio 
remnant within the Vela supernova remnant. The fact that the pulsar nebula is displaced from 
the pulsar can be explained in our model without the need for directed power from the pulsar. 
Asymmetries in the placement of radio PWNae relative to their pulsars are likely to be common in 
this phase because the interstellar medium is generally inhomogeneous. 

Support for this work was provided in part by NASA grants NAG5-8130 and NAG-7153. The 
hydrodynamic simulations were computed on an IBM SP at the North Carolina Supercomputing 
Center. 
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